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ABSTRACT 

This proposed two-year research project was to involve development of an analytical model, a 
numerical algorithm for its integration, and a software for the analysis of a solidification process 
under the influence of electric and magnetic fields in microgravity. Due to the complexity of the 
analytical model that was developed and its boundary conditions, only a preliminary version of 
the numerical algorithm was developed while the development of the software package was not 
completed. 


SUMMARY OF TECHNICAL ACCOMPLISHMENTS 

We have developed and refined a new analytical formulation of combined electro-magneto- 
hydrodynamics (EMHD) and have shown the inconsistencies and shortcomings of the existing 
separate electro-hydrodynamic (EHD) and magneto-hydrodynamic (MHD) theories. The new 
model consists of a system of at least twelve coupled non-linear partial differential equations in 
case of a three-dimensional flow. The model is valid for multi-dimensional, unsteady, viscous 
fluid flows involving electrically charged particles and electric linear polarization and 
magnetization effects. All interaction between externally applied and internally induced electric 


2 


and magnetic fields have been incorporated in the model. We have also developed a fully 
conservative vector operator form of the EMHD system that is suitable for a direct discretization 
and numerical integration. Derivations of characteristic and non-reflecting open boundary 
conditions were also performed, thus completing the analytical part of the research effort. A 
numerical algorithm for iterative integration of the discretized fully conservative EMHD system 
in boundary-conforming non-orthogonal curvilinear coordinate system was incompletely 
developed that is based on dual time-stepping. Development of an accompanying computer code 
for the analysis of unsteady three-dimensional EMHD flows is in its initial debugging stages. 

SUMMARY OF MEASURABLE ACCOMPLISHMENTS 

INDUSTRY SUPPORT: one 

ALCOA Research Center 

• provided a Faculty Research Fellowship to the P. I. to work on inverse problems 

GOVERNMENT SUPPORT: two 

NASA Ames Research Center NAS facility 

• provided free computing time on a Cray C-90 computer 
NASA Marshall Space Flight Center 

• provided a three day training for a graduate student working on an experiment 

NUMBER OF STUDENTS FUNDED: one 

Sean R. Lynn: "A Unified Theory and Simulation of Unsteady Electro-Magneto- 

Hydrodynamics," M.Sc. degree candidate. 

(Student accepted a permanent employment offer from Allison Gas Turbines 
Company on October 1, 1996. Does not intend to return for completion of his 
degree.) 

NUMBER OF RELATED B.Sc. HONORS THESIS COMPLETED: two 

Bates, Craig: "Inverse Determination of Locations and Strengths of Electric Impulses 
Inside a Human Heart Based on Chest Surface Measurements of Electric 
Potential", B. Sc. Honors Thesis, Dept, of Eng. Science & Mechanics, The 
Pennsylvania State University, May 1996. 

Gross, Chris: "Feasibility Study of a Magnetohydrodynamic Blood Pump", B. Sc. 

Honors Thesis, Dept, of Eng. Science & Mechanics, The Pennsylvania State 
University, May 1996. 

NUMBER OF RELATED M.Sc. THESIS COMPLETED: one 

Craig Bates: "Forward and Inverse Electro-Cardiographic Calculations on a Multidipole 
Model of Human Cardiac Electrophysiology”, Dept, of Eng. Science & 
Mechanics, The Pennsylvania State University, August 1997. 

NUMBER OF RELATED B.Sc. HONORS THESIS IN PROGRESS: one 

James Tao: “Inverse Determination of Magnet Strengths Generating a Desired Magnetic 
Field Pattern”, Dept, of Eng. Science & Mechanics, The Pennsylvania State 
University, expected May 1998. 
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NUMBER OF RELATED PAPERS PUBLISHED IN REFEREED JOURNALS: five 

Dulikravich, G. S. and Lynn, S. R., “Unified Electro-Magneto-Fluid Dynamics (EMFD): 
A Survey of Mathematical Models,” International Journal of Non-Linear 
Mechanics, Vol. 32, No. 5, September 1997, pp. 923-932. 

Dulikravich, G. S. and Lynn, S. R., “Unified Electro-Magneto-Fluid Dynamics (EMFD): 
Introductory Concepts,” International Journal of Non-Linear Mechanics, Vol. 
32, No. 5, September 1997, pp. 913-922. 

Martin, T. J. and Dulikravich, G. S., “Inverse Determination of Boundary Conditions in 
Steady Heat Conduction with Heat Generation,” ASME Journal of Heat 
Transfer, Vol. 1 18, No. 3, August 1996, pp. 546-554. 

Choi, K. -Y. and Dulikravich, G. S., “Acceleration of Iterative Algorithms on Highly 
Clustered Grids,” AIAA Journal, Vol. 34, No. 4, April 1996, pp. 691-699. 

Martin, T. J. and Dulikravich, G. S., “Finding Unknown Surface Temperatures and Heat 
Fluxes in Steady Heat Conduction,” IEEE Transactions on Components, 
Packaging and Manufacturing Technology (CPMT) - Part A, Vol. 18, No. 3, 
September 1995, pp. 540-545. 

NUMBER OF RELATED JOURNAL ARTICLES IN PRESS: two 

Martin, T. J. and Dulikravich, G. S., “Inverse Determination of Steady Convective Local 
Heat Transfer Coefficients,” ASME Journal of Heat Transfer. 

Martin, T. J. and Dulikravich, G. S., “Non-Iterative Determination of Temperature- 
Dependent Heat Conductivities,” ASME Journal of Heat Transfer. 

NUMBER OF RELATED BOOK CHAPTERS AND PROCEEDINGS: four 

Dulikravich, G. S.: “Design and Optimization Tools Development”, Chapters no. 10- 
15 in New Design Concepts for High Speed Air Transport, (editor: H. 
Sobieczky), Springer, Wien/New York, 1997, pp. 159-236. 

Dulikravich, G. S. and Martin, T. J.: "Inverse Shape and Boundary Condition 

Problems and Optimization in Heat Conduction", Chapter no. 10 in Advances 
in Numerical Heat Transfer - Volume 1 (editors: W. J. Minkowycz and E. M. 
Sparrow), Taylor and Francis, November 1996, pp. 381-426. 

Proceedings of Symposium on Inverse Design Problems in Heat Transfer and Fluid 
Flow, Editors: G. S. Dulikravich and K. A. Woodbury, ASME National Heat 
Transfer Conference, Baltimore, MD, August 10-12, 1997, ASME HTD-Vol. 340, 
Volume 2. 

Proceedings of Symposium on Developments in Electrorheological Flows-1995, 

Editors: D. A. Siginer and G. S. Dulikravich, ASME WAM'95, San Francisco, 
CA, November 12-17, 1995, ASME FED- Vol. 235, MD-Vol. 71. 

NUMBER OF RELATED BOOK CHAPTERS AND PROCEEDINGS IN PRESS: one 

Dulikravich, G. S., “Electro-Magneto-Hydrodynamics and Solidification,” Chapter 
no. 9 in Advances in Non-Newtonian Flows and Rheology, (editors: D. A. 
Siginer, D. De Kee and R. P. Chhabra), Elsevier Publishers, Spring 1998. 

REFEREED JOURNAL ARTICLES IN PREPARATION: one 

Dulikravich, G. S., "Unified Electro-Magneto-Hydrodynamics: A Mathematical Model 
and Boundary Conditions", Physics of Fluids. 
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PAPERS PRESENTED AT TECHNICAL MEETINGS: twelve 

Dulikravich, G. S. and Jing, Y. -H., “Boundary Conditions for Electro-Magneto- 

Hydrodynamics,” Symposium on Rheology and Fluid Mechanics of Non- 
Linear Materials, Editors: S. G. Advani and D. A. Siginer, ASME IMECE'97, 
Dallas, TX, Nov. 16-21, 1997, ASME FED-Vol. 243/MD-Vol. 78, pp. 101-117. 

Dulikravich, G. S. and Martin, T. J., “Inverse Determination of Steady Local Convective 
Heat Transfer Coefficients,” Symposium on Inverse Design Problems in Heat 
Transfer and Fluid Flow, Editors: G. S. Dulikravich and K. A. Woodbury, 
ASME National Heat Transfer Conference, Baltimore, MD, August 10-12, 1997, 
ASME HTD-Vol. 340, Volume 2, pp. 1 5 1 - 1 58. 

Dulikravich, G. S. and Martin, T. J., “Non-Iterative Inverse Determination of 
Temperature-Dependent Heat Conductivities,” Symposium on Inverse Design 
Problems in Heat Transfer and Fluid Flow, Editors: G. S. Dulikravich and K. 
A. Woodbury, ASME National Heat Transfer Conference, Baltimore, MD, 
August 10-12, 1997, ASME HTD-Vol. 340, Volume 2, pp. 141-150. 

Dulikravich, G. S. and Martin, T. J., “Inverse Determination of Boundary Conditions in 
Field Problems, (Invited Lecture), Advanced Technology in Experimental 
Mechanics - ATEM97, Editor: Y. Morimoto, Wakayama City, Osaka, Japan, 
July 25-26, 1997, pp. 83-88. 

Martin, T. J. and Dulikravich, G. S., “Inverse Determination of Boundary Conditions in 
Multidomain Heat Conduction Problems,” (Invited Lecture), BETECH '97 - 9th 
International Conference on Boundary Element Technology, Editor: J. 
Frankel, Knoxville, TN, April 9-11, 1997, pp. 99-1 10. 

Dulikravich, G. S. and Martin, T. J., “Inverse Determination of Temperatures and Heat 
Fluxes on Surfaces of 3-D Objects,” Pan-American Congress of Applied 
Mechanics (PACAM-V), San Juan, Puerto Rico, January 2-4, 1997, in Applied 
Mechanics in the Americas, Editors: M. Rysz, L. A. Godoy and L. E. Suarez, 
Vol. 5, pp. 133-136. 

Dulikravich, G. S. and Jing, Y. -H., “Fully Conservative Forms of a System of Unified 
Electro-Magneto-Hydrodynamic Equations,” Symposium on Rheology and 
Fluid Mechanics of Nonlinear Materials, Editors: D. A. Siginer and S. G. 
Advani, ASME International Mechanical Engineering Congress and Exposition, 
Atlanta, GA, November 17-22, 1996, AMD- Vol. 217, pp. 309-3 14. 

Dulikravich, G. S., "An Unified Theory of Electro-Magneto-Gasdynamics", Hypersonic 
MAGLEV Group Meeting, National High Magnetic Field Lab., Tallahassee, FL, 
December 12, 1995. 

Dulikravich, G. S. and Lynn, S. R., "Electro-Magneto-Hydrodynamics: Part 2 - A Survey 
of Mathematical Models", Symposium on Developments in Electrorheological 
Flows-1995, Editors: D. A. Siginer and G. S. Dulikravich, ASME WAM'95, San 
Francisco, CA, November 12-17, 1995, ASME FED-Vol. 235, MD-Vol. 71, pp. 
59-70. 

Dulikravich, G. S. and Lynn, S. R., "Electro-Magneto-Hydrodynamics: Part 1 - 
Introductory Concepts", Symposium on Developments in Electrorheological 
Flows-1995, Editors: D. A. Siginer and G. S. Dulikravich, ASME WAM'95, San 
Francisco, CA, Nov. 12-17, 1995, ASME FED-Vol. 235, MD-Vol. 71, pp. 49-58. 
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Dulikravich, G. S., Choi, K. -Y. and Lee, S., “Magnetic Field Control of Vorticity in 
Steady Incompressible Laminar Flows,” Symposium on Developments in 
Electrorheological Flows and Measurement Uncertainty 1994, ASME 
WAM'94, Editors: D. A. Siginer, J. H. Kim, S. A. Sheriff and H. W. Colleman, 
Chicago, IL, November 6-11,1994, ASME FED-Vol. 205/AMD- Vol. 190, pp. 
125-142. 

Dulikravich, G. S. and Martin, T. J., “Inverse Problems and Design in Heat Conduction,” 
2nd IUTAM International Symposium on Inverse Problems in Engineering 
Mechanics, Editors: H. D. Bui, M. Tanaka, M. Bonnet, H. Maigre, E. Luzzato 
and M. Reynier, Paris, France, November 2-4, 1994, A. A. Balkema, Rotterdam, 
1994, pp. 13-20. 

OTHER PERTINENT TECHNICAL PRESENTATIONS: twenty 

Invited Lecture, Graduate School of Eng., Kyoto University, Kyoto, Japan, August 1995.: 
"Inverse Problems and Optimization in Heat Transfer" 

Invited Lecture, Mechanical Eng. Dept., Shinshu University, Nagano, Japan, August 1995.: 
"Inverse Problems and Optimization in Heat Transfer" 

Invited Lecture, Ebara Research Company, Ebara Company, Kanagawa, Japan, August 1995.: 
"Inverse Problems and Optimization in Fluid Mechanics" 

Invited Lecture, National Aerospace Laboratory - NAL, Tokyo, Japan, August 1995.: 
"Multidisciplinary Inverse Problems and Optimization" 

Invited Lecture, Mechanical Eng. Dept., Ashikaga Inst, of Tech, Ashikaga, Japan, Aug 1995.: 
"Electro-Magneto-Hydrodynamics and Solidification" 

Invited Lecture, Mechanical Eng. Dept., University of Tokyo, Tokyo, Japan, July 1995.: "Inverse 
Problems and Optimization in Fluid Mechanics" 

Invited Lecture, Ishikawajima-Harima Heavy Industries R & D, Tokyo, Japan, July 1995.: 
"Multidisciplinary Inverse Problems and Optimization" 

Invited Lecture, Fundamental Research Labs, NEC Corp., Tsukuba, Japan, July 1995.: "Electro- 
Magneto-Hydrodynamics and Solidification" 

Invited Lecture, Dept, of Aero. & Space Eng., Tohoku University, Japan, July 1995.: "Inverse 
Problems and Optimization in Fluid Mechanics" 

Invited Lecture, Toshiba Corp. R & D Center, Kawasaki, Japan, July 1995.: "Multidisciplinary 
Inverse Problems and Optimization" 

Invited Lecture, Mechanical Eng. Lab., Hitachi, Ltd., Tsuchiura, Japan, July 1995.: 
"Multidisciplinary Inverse Problems and Optimization" 

Invited Lecture, Mechanical Faculty, Nat. Tech. Univ. of Athens, Athens, Greece, June 1995.: 
"Multidisciplinary Inverse Problems and Optimization" 

Invited Lecture, Mechanical Faculty, Aero. Inst., Univ. of Belgrade, Yugoslavia, May 1995.: 
"Inverse Problems and Optimization in Fluid Mechanics" 

Lecture, Aerospace Eng. Dept., Pennsylvania State Univ., University Park, PA, March 1995.: 
"Inverse Problems and Optimization in Heat Transfer" 

Lecture, Center for Theor. and Comput. Materials Science, NIST, Gaithersburg, Feb. 1995.: 
"Electro-Magneto-Hydrodynamics and Solidification" 

Invited Lecture, Mechanical Eng. Dept., Univ. of Minnesota, Minneapolis, MN Jan. 1995.: 
"Inverse Problems and Optimization in Heat Transfer" 

Invited Lecture, Mechanical Eng. Dept., University of Pittsburgh, Pittsburgh. PA, Jan. 1995.: 
"Electro-Magneto-Hydrodynamics and Solidification" 

Invited Lecture, ALCOA Technical Center, ALCOA Center, PA, August 1996.: "Electro- 
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Magneto-Hydrodynamics and Solidification" 

Invited Lecture, Mech. Eng. Dept., California State University, Fullerton, CA, July 1996.: 

"Inverse Problems and Optimization in Heat Transfer" 

Invited Lecture, Mech. Eng. Dept, The Johns Hopkins Univ., Baltimore, MD, Feb. 1996.: 

"Inverse Problems and Optimization in Heat Transfer" 

PERTINENT PROFESSIONAL MEETINGS AND SESSIONS ORGANIZED: sixteen 

“Forum on Functional Fluids”, forum co-organizer at 1999 Joint ASME/JSME Fluids 
Engineering Conference, San Francisco, CA, July 18-23, 1999. 

“Multidisciplinary Inverse Problems and Optimization in Heat Transfer”, co-organizer 
and co-chairman of the Symposium at ASME IMECE98, Anaheim, CA, 
November 15-20, 1998. 

“Elastic Fluids,” co-chairman of the session at ASME IMECE, Dallas, TX, November 
16-21, 1997. 

"Inverse Design Problems in Heat Transfer and Fluid Flow,” co-organizer and co- 
chairman of the Symposium at ASME National Heat Transfer Conference, 
Baltimore, MD, August 10-12, 1997. 

"Advanced Technology in Experimental Mechanics- ATEM97,” invited speaker, session 
chair, and member of the International Program Committee, Wakayama City, 
Osaka, Japan, July 25-26, 1997. 

"BETECH '97 - 9th International Conference on Boundary Element Technology," 
member of the Scientific Advisory Committee, Knoxville, TN, April 9-11, 1997. 

"Heat Transfer", session chairman at the Pan-American Congress of Applied Mechanics 
(PACAM-V), San Juan, Puerto Rico, January 2-4, 1997. 

"Rheology and Fluid Mechanics of Nonlinear Materials IV: Complex Flows", session 
vice-chairman at ASME IMECE'96, Atlanta, GA, Nov. 17-22, 1996. 

"Second International Conference on Inverse Problems in Engineering: Theory and 
Practice", member of the Scientific Advisory Committee, Nantes, France, June 
1996. 

"BETECH '96 - 9th International Conference on Boundary Element Technology," 
member of the Scientific Advisory Committee, Maui, Hawaii, April 24-26, 1996. 

"3rd International Symposium on Magnetic Suspension Technology", session chairman, 
Tallahassee, FL, December 13-15, 1995. 

"Symposium on Electrorheological Flows - IH", session co-organizer, ASME WAM'95, 
San Francisco, CA, November 12-17, 1995. 

"The Seventh Inverse Problems in Engineering Seminar", member of the Organizing 
Committee, Columbus, OH, June 12-13, 1995. 

"PACAM IV- Pan-American Congress of Applied Mechanics," member of the 
Organizing Committee, Buenos Aires, Argentina, January 3-6, 1995. 

"Symposium on Inverse Problems in Mechanics - III", session co-chairperson, ASME 
WAM'94, Chicago, IL, November 6-11, 1994. 

"Symposium on Inverse Problems in Engineering Mechanics ISIP'94," member of the 
International Scientific Committee, November 2-4, 1994, Paris, France. 

HONORS AND AWARDS: two 

Fellow, American Society of Mechanical Engineers (ASME), Nov. 1996. 

ALCOA Foundation Faculty Fellow Research Award (July 1,1996 - June 30, 1998) 


ELECTRO-MAGNETO-HYDRODYNAMICS 
AND SOLIDIFICATION 

G. S. Dulikravich 

Aerospace Engineering Department, The Pennsylvania State University, 
University Park, Pennsylvania 16802, USA 


1. INTRODUCTION 

Fluid flow influenced by electric and magnetic fields has classically been 
divided into two separate fields of study: electro-hydrodynamics (EHD) 

studying fluid flows containing electric charges under the influence of an 
electric field and no magnetic field, and magneto-hydrodynamics (MHD) 
studying fluid flows containing no free electric charges under the influence of a 
magnetic field and no electric field. Traditionally, this division was necessary 
to reduce the extreme complexity of the coupled system of Navier-Stokes, 
Maxwell’s and constitutive equations describing combined electro-magneto- 
hydrodynamic flows. Recent advances in numerical techniques and computing 
technology, as well as fully rigorous theoretical treatments, have made analysis 
of combined electro-magneto-hydrodynamic flows well within reach. A survey 
of electro-magnetics and the theory describing combined electro-magneto- 
hydrodynamic (EMHD) flows is presented with an emphasis on describing the 
intricacies of the mathematical models and the corresponding boundary 
conditions for fluid flows involving polarization and magnetization. This 
survey concludes with a presentation of EHD and MHD flow models involving 
solidification. ° 


NOMENCLATURE 

b = electric charge mobility coefficient, kg A s 2 

B = magnetic flux density vector, kg A* 1 s' 2 
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— average rate of deformation tensor, s’ 1 

= electric charge diffusion coefficient, m 2 s' 1 
= electric displacement field vector, A s m' 2 
= total energy per unit mass, m 2 s' 2 

= electric field vector, kg m s' 3 A' 1 , or V m' 1 
= electromotive intensity vector, kg m s' 3 A' 1 
= mechanical body force vector per unit mass, m s' 2 
= acceleration due to gravity, m s' 2 
= heat source or sink per unit mass, m 2 s' 3 
= magnetic field intensity vector, A m' 1 
= electric current density vector, A m' 2 
= electric conduction current vector, A m' 2 
= electric drift current vector, A m' 2 
= total magnetization vector per unit volume, A m' 1 
= magnetomotive intensity vector per unit volume, A m 1 
= pressure, kg m' 1 s' 2 

= total polarization vector per unit volume, A s m' 2 
= total or free electric charge per unit volume, A s m' 3 
= heat flux vector, kg s' 3 

= entropy per unit mass, m 2 kg' 1 K ‘ s' 2 
= absolute temperature, K 
= fluid velocity vector, m s' 1 


GREEK SYMBOLS 


a 

P 


= volumetric thermal expansion coefficient, K' 1 
= Chorin’s (1967) artificial compressibility coefficient 


£ =dielectric constant (electric permittivity), kg' 1 m' 3 s 4 A 

e 0 = 8.854 x 10 12 = vacuum electric permittivity, kg' 1 m' 3 s 4 A 2 

e r = e / e 0 = relative electric permittivity 

K = thermal conductivity coefficient, kg m s' 3 K' 1 

a = electric conductivity coefficient, kg' 1 m' 3 s 3 A 2 

p = fluid density, kg m' 3 

1 V = + M-v 2 |(^ • v) = Newtonian viscous stress tensor, kg m' 1 s' 2 

l = electromagnetic stress tensor, kg m' 1 s' 2 



Ho =4rcx 10 " 7 

M-r =H/Ho 
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= stress tensor (viscous plus electromagnetic), kg m' 1 s’ 2 
= magnetic permeability coefficient, kg m A' 2 s' 2 

= magnetic permeability of vacuum, kg m A' 2 s' 2 
= relative magnetic permeability 
= shear coefficient of viscosity, kg m' 1 s' 1 
= second coefficient of viscosity, kg nT 1 s' 1 
= electric susceptibility 
= magnetic susceptibility 
= electric potential, V 
= viscous dissipation function, kg m 's 3 


2 . BACKGROUND 

The scientific field of study that analyzes the ability of electro-magnetic 
fields to influence fluid flow-field and heat transfer has been investigated for 
decades. The equations that are most often used to model this phenomena 
consist of the system of Navier-Stokes equations for fluid motion coupled with 
Maxwell’s equations of electro-magnetics augmented with the material 
constitutive relations. The field studying these flows is often called electro- 
magneto-dynamics of fluids [1], electro-magneto-fluid dynamics (EMFD) [2-5], 
electro-magneto-hydrodynamics [6], magneto-gas-dynamics and plasma 
dynamics [7], or the electro-dynamics of continua [8-10]. The full system of 
governing equations has, until recently, been far too difficult to solve because 
Navier-Stokes system becomes very complex when modeling flows involving 
turbulence, chemical reactions, multiple phases, non-Newtonian effects, etc. 
When coupled with Maxwell’s equations, the complexity of the combined 
EMHD system is raised by orders of magnitude. To reduce this complexity, the 
analytical modeling has traditionally been divided [11] into flows influenced 
only by externally applied electric fields acting upon electrically charged 
particles in the fluid, and flows influenced only by externally applied magnetic 
fields without electric charges in the fluid. The former are called Electro- 
Hydrodynamic (EHD) flows [12] and the latter Magneto- Hydrodynamic (MHD) 
flows [13]. More recently, rigorous continuum mechanics treatments of EHD 
[14] and unified EMHD flows [9,10] have been developed. These continuum 
mechanics approaches are limited to non-relativistic, quasi-static or relatively 
low frequency phenomenon [15-17]. 



This chapter should provide an introductory survey of the background theory 
to allow implementation of numerical analysis of unified EMHD flows and of 
classical MHD and EHD flows with addition of liquid/solid phase change. An 
overview of electro-magnetic theory with concentrated effort placed on 
descriptions of the electric and magnetic fields and electric charges and currents 
will be made to provide a physical understanding of the field-material 
interactions causing polarization and magnetization effects. The system of 
equations governing the unified EMHD theory and the corresponding boundary 
conditions will be presented together with its fully conservative form that is 
ready for numerical discretization. 


3. POLARIZATION AND GAUSS’ LAW 

Charge polarization is created when electric charges of opposite signs are 
separated by a distance. Although many references define several sources of 
polarization [18], there are essentially two main sources of polarization: natural 
and induced [13]. Natural polarization arises from natural dipoles and charged 
particles. An example of a natural dipole is a water molecule which has a 
geometry such that the centers of positive charges and negative charges do not 
coincide. Since the molecules are allowed to move freely and orient randomly, 
water will not have polarization on a continuum level. Now consider the fluid 
water as it is frozen with an applied electric field. An induced polarization will 
be created by the electric field by inducing an initial charge separation in 
neutral particles [19], by causing greater charge separation within the 
molecules, and by causing molecular alignment with the applied electric field in 
case of natural dipoles [19]. Once locked in the ice crystal structure, the water 
molecules will no longer be able to change their position or orientation. 
Consequently, even after the electric field is removed, the ice will still have 
polarization on a continuum level since the polarization caused by the electric 
field aligning the water molecules was literally frozen into the ice. 

From this example it may seem that there is no reason, when dealing with 
fluids, to consider natural polarization. This, however, would be an erroneous 
assumption. Though the natural polarization may show no continuum effects 
without the presence of an electric field, in an electric field the total 
polarization, P, combines both the induced polarization due to the electric field 
and the natural polarization of the molecules which are now aligned by the 
electric field [13, p.22]. 



If polarization is assumed to be a linear function of the steady or relatively 
low frequency electric field, then it can be defined as 

P = e 0 X E (E + yxB)=e p (E + yxB)=8pE (1) 

The electric displacement vector then becomes [19, p.164] [9, p.178]. 

D = e 0 E + P = e 0 (l + % E )E + £o% E yxB = s 0 e r E + e p yxB = eE + e p vxB (2) 

where the material property, % E , is the dielectric susceptibility. It is typically 
obtained experimentally [20, p.86] and could be a function of frequency. 

Electric charges come in two types: free and bound. Free charges arise from 
electrons in the outer or free atomic shells and from ions. Bound charges are 
those arising from the molecular geometry and displacement of atomic inner 
electron shells [13, p.21]. Gauss’s law for a linearly polarizable medium then 
becomes [13, p. 22] 


V-D = q 0 
or 


(3) 


^•(e 0 E + E)— V(eE + s p yxB)=q 0 ^ 

At this point it is important to note that q 0 multiplied with the charged particle 
drift velocity, y d , creates the convection or drift electric current, J d [13, p.67], 

while polarization current, Ip, is defined as the variation of the total 
polarization with respect to time [19, p.121 and p.147]. 


4. MAGNETIZATION AND AMPERE-MAXWELL’S LAW 

If the material in question may be considered linear, that is, if the 
magnetization is a function of one material property and the strength and 
direction of the applied magnetic field, then the magnetization is defined as \9 
P-178] [19, p. 164] [20, p.92-96] [21, p.371-377] 


M = M + vxP = — -h r B 

H, & + *“)- 


( 5 ) 



In addition to the electric currents arising from magnetization and direct charge 
motion, other phenomenological currents have been observed and must be taken 
into account when defining the total current, J [9, p.162-1 63], Introducing the 
effects of magnetization and polarization and rearranging constants, the 
Ampere-Maxwell’s law of electrodynamics may be rewritten as [13, p.30] 


VxB = n 0 (vxM+J d +I„+^I'| 

l dt , 


( 6 ) 


Magnetization and magnetic field vectors are often combined to form the 
magnetic field strength vector, H, defined as 


H = — -M 
Ho 


(7) 


The total current, J, is defined as the sum of the apparent magnetization current, 
VxM, charge drift current, J d , and phenomenological polarization currents, 
sip [13, p.26] since the contribution to the magnetization current by intrinsic 

magnetization is zero. The Ampere-Maxwell’s law for polarizable, magetizable 
media can therefore be written as [19, p. 132] 

f- Vx H = -I (8) 

Detailed descriptions of these equations can be found in any number of texts 
[19, 20,21]. 


5. A MODEL OF UNIFIED ELECTRO-MAGNETO-GASDYNAMICS 
(EMGD) 

The full system of equations governing unified EMGD flows consists of the 
Maxwell’s equations governing electro-magnetism, the Navier-Stokes equations 
governing compressible fluid flow, and constitutive equations describing 
material behavior. Assuming a single-phase fluid and only one type of charged 
particles in the fluid, this set has a minimum of 12 partial differential equations 
that contains 13 unknowns: p, q o , T, p, and the three vector components of v, E, 
and B, respectively. The thirteenth equation is the equation of state for the 


fluid. The foundations of the electro-magneto-gasdynamic (EMGD) theory 
were formulated by Eringen and Maugin [9,10] and are based on continuum 
mechanics [22-25]. The rigor with which the constitutive, force, and energy 
terms were derived leads to a model more complete and robust than any of those 
found in classical literature [8,7,1, 1 1-13,1 8-21]. 

Dulikravich and Jing [6,26] have shown that a compact vector form of the 
unified EMGD system can be written as a combination of the Maxwell’s 
electro-magnetic subsystem and the Navier-Stokes fluid flow subsystem. 

The Maxwell’s subsystem (consisting of seven PDE’s) is composed of 
Ampere-Maxwell’s law for polarizable and magnetizable medium 

dD V7 TT 

-^-VxH = -I (9 ) 


which can also be written as 
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Faraday’s law 


— +VxE=0 

at 


( 10 ) 


(ii) 


and conservation of electric charges 

^tv-l =0 

at 

which is a combination of Gauss’ law 


( 12 ) 


V D = q 0 (13) 

and the Ampere-Maxwell’s law. Conservation of magnetic flux 

VB = 0 (14) 


is also a part of the Maxwell’s subsystem, but is not solved for explicitly. 


The second part of the unified EMGD is the viscous, compressible flow 
Navier-Stokes subsystem consisting of five PDE’s and an equation of state of a 
perfect gas. It is composed of conservation of mass equation 

^+V-(py) = 0 (15) 

and a conservation of linear momentum (including electromagnetic effects) 

+ V • (ypv - x) - V • (v(P x B )) - V • ((B ■ M)I + (E ■ P)|) = S v (16) 


Here, I is the identity (unity) tensor and SMs a vector of source terms. The 
following dyadic identities were used in equation (16) 

(VB)M = V((BM|)-(VM)B (17) 

(VE)P = V((EP|)-(VP).E (18) 

Conservation of energy equation is also a part of the Navier-Stokes subsystem 
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It can be replaced by the Clausius-Duheim entropy inequality [9,2,4] 
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The viscous stress tensor for a non-linear fluid is given as 
l =2p v d + |j. v2 I(V-y) + a 1 d 2 


( 21 ) 



In the case of a media with non-linear physical properties, the unified EMGD 
formulation for the electric conduction current and the heat flux can be 
expressed as [9, p. 161- 1 62], 

I c = o,E+o 2 d| + o 3 d 2 E + o 4 VT + a 5 d- VT + a 6 d 2 • VT + g 7 ExB 

+ a 8 (d-(ExB)-(d-E)xB) + a 9 VTxB (22) 

+ ^io(d(VTxB)-(dVT)xB) + a n (B-E)B + a 12 (BVT)B 

q = K,E + K 2 d-E + K 3 d 2 E + K 4 VT + K 5 d-VT + K 6 d 2 -VT + K 7 ExB 
+ k 8 (d • (E x B) - (d • E) x B) + k 9 VT x B (23) 

+ K i 0 (d ■ ( VT x B) - (d • VT) x B) + k, , (B • E)B + k 12 (B • VT)B 

The electro-magneto-thermal stress tensor for a non-linear fluid is [9, p.177- 
178] 

x EM =a 2 E®I+a 3 B®B + a 4 VT®VT+a 5 (E(8>dE) s 
+ a 6 (E®d 2 i) s +a 7 (VT<8)d- VT) $ + otg(vT<8>d 2 ■ VT^ 

+ a 9 (dW-W-d)+a 10 WdW + a 11 ^ 2 -W-Wd 2 ) 

+ a 12 fe-d-W 2 -W 2 -d-w)+a 13 (E(8)VT) s +a 14 (W-E(g)EW) s 

+ a 15 (E® WE) s +a 16 (w-E® W 2 -E^ +a 17 (W-(E(8) VT- VT(8 >E)) s 
+ a 18 d • (E ® VT - VT (8) E)- a 18 (E ® VT - VT ® E)- d 

where W = = £y k B k , while the subscript s indicates symmetrization. 

Expressions for total polarization, P, and magnetization, M, of non-linear media 
can be modeled with expressions of similar complexity [9, p. 175]. 

In these formulas, a., Gj and k are the physical properties of the media. Most 
of these coefficients are still unknown although their exploitation can offer 
potentially significant benefits in applications involving interacting electric, 
magnetic, thermal, and stress fields. This theory is valid for the frequencies of 
the electric and the magnetic fields that are less than approximately 1 kHz and 
for fluid speeds considerably less than the speed of light [14-17]. For higher 
frequencies, certain physical properties become functions of the frequencies. 
For higher speeds, relativistic effects will have to be taken into account. 


6. CONSERVATIVE FORMS OF ELECTRO-MAGNETO- 
HYDRODYNAMIC (EMHD) SYSTEM 

A necessary condition that an iterative numerical solution of the EMGD 
system will converge to the exact solution of the analytical EMGD system as 
the computational grid is infinitely refined, requires that the EMGD system 
must be rewritten in a fully conservative (divergence-free) form. This is 
especially needed if strong gradients of dependent variables are expected to 
exist in the solution domain. The fully conservative forms can then be used 
directly in the finite difference, finite volume, or finite element discretization of 
the EMGD system and its iterative integration process. 

In the following derivations, it will be assumed that the fluid is 
incompressible, homocompositional, that it has linear polarization and linear 
magnetization properties, and that the frequencies of the applied electric and 
magnetic fields are less than approximately 1000 Hz for this mathematical 
model to be realistic. These are the only assumptions to be used in this model 
which will be referred to as a unified electro-magneto-hydrodynamics (EMHD). 

A fully conservative EMHD system in a vector operator form is given as [6] 
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For simplicity of notation we can define the following terms as [6,26] 
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If we now assume that the fluid which is subjected to applied electric and 
magnetic fields is of Newtonian type and if we allow only for linear 
polarization (equation 1) and linear magnetization (equation 5), the constitutive 
relations for the electric conduction current and for the heat flux vector become 
[9, p. 173- 174] 

J c =a 1 (E+yxB) + a 4 VT+a 7 (E + yxB)xB 

+ G 9 VTxB + a n (B-(E + yxB))B + a 12 (B-VT)B 


(38) 


( 39 ) 


q = Kj(E + vxB) + k 4 VT+k 7 (E + vxB)xB 
+ k 9 VT x B + k, j (B • (E + v x B))B + k 12 (B • VT)B 

Then, the EMHD source terms can be given in a compact vector form [6,26] as 

S E =--(!+£) 

e o 

S v =f + ^[q 0 E + (VxE)xP-(VM)B-(VP>E + (l + P 1 )xB] 

S'=h + i(E + vxB)[(vV^ + J c +p i ]-iir x M B((vV)B-VxE) 

Notice that these source terms have been formulated in such a way as not to 
contain explicit time derivatives [6,26]. 
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6.1 Fully conservative Cartesian form of the EMHD system 

The EMHD system of equations (equations 25-30) can now be written in a general 
conservative form in terms of (x,y,z) orthogonal coordinate system as 
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Here, the solution vector of unknown quantities is given as 








(44) 


where the asterisk symbol designates transpose of a vector. The vector of 
source terms (those terms that do not contain divergence operator) is given as 




0, 0, 0, 0, S x v , Si, S 





(45) 


In equation (44), Chorin’s [27] artificial compressibility coefficient, (3, was used 
to create the unsteady term in the mass conservation since physical unsteady 
term does not exist in the mass conservation for incompressible fluids. 

By combining equations (5), (7), (1), and (31), the Cartesian components of 
the magnetic field intensity vector can be defined [6,26] as 

H * =F B x +e p v y (E z + v x B y - v y B x )-e p v z (E y + v z B x - v x B z ) 

Hy = F B y +e p v z (E x + v y B z -v z B y )-e p v x (E z +v x B y - v y B x ) (46) 

H z = p B z +e P v x (E y + v z B x - v x B z )-8 p v y (E x +v y B z -v z B y ) 

The flux vectors in equation (43) can then be defined as 
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Here, we have written components of (PxB) as 


Nj B =PyB z -P z B y 

Ny B =P z B x -P x B z 

N BB =P x B y -P y B x 


In addition, we have defined the terms 
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Components of the electric current vector, J, were defined as 
Jx=v x q 0 +^ L P x +a 4 ^+^-N p x B +a 9 (^B z -|tB v ) 
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and heat flux vector components were defined as 


is lf ltI( E4N,' + K,eB I ~» 

dx Ep dy dz 


aT " }r B,)t^N 1 A +% N 11 B, 


+ | + | N ” + K ’€ B « -f B 2 ) + ^N„B y + K|2 N„B, (56) 


q 2 =^-P I + K 4 £+|lN™+K 9 (^By— ^BJ+^-N b 1 ,B 2 + K 12 N„B 


dx y dy 



7. CHARACTERISTIC-BASED INFLOW AND OUTFLOW 
BOUNDARY CONDITIONS 

For most boundary value problems of electro-magneto dynamics, jump 
conditions are exclusively used [9,28] to formulate solid wall boundary 
conditions where a discontinuity occurs. At the inflow and outflow boundaries 
where no surface or line discontinuities exist, an alternative approach based on 
conservation law for continuous surfaces or lines become necessary. 
Characteristic boundary condition formulation [29,30], which starts from a 
characteristic form of the EMHD system, will be sketched here since it leads to 
non-reflecting boundary condition formulation [31-36,26]. To find the 
characteristic boundary conditions, it is first necessary to determine analytical 
expressions for all eigenvalues of the characteristic system. The most common 
approach is to use one of the symbolic programming languages software (LISP, 
MACSIMA) in order to determine analytical expressions for each eigenvalue. 
Since these software packages cannot be used for systems that have more than 
five coupled partial differential equations, in the case of a complete EMHD 
system which has twelve coupled partial differential equations, it is impossible 
to find the eigenvalues using available symbolic programming software. 

Consequently, we will use an alternative approach in which we will divide 
the unified EMHD system into a Maxwell’s subsystem and the Navier-Stokes 
subsystem [33]. Each of these two subsystems will then be analyzed separately 
by finding the analytical expressions for its eigenvalues by hand. 

7.1 Characteristic-based boundary conditions for Maxwell’s subsystem 

For example, characteristic treatment of the Maxwell’s subsystem can be 
formulated by rewriting the fully conservative Maxwell’s subsystem 

^Qem , dE EM , dF E M , d(jEM _ g 
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in a non-conservative (characteristic) form as 
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In order to perform characteristic analysis for Maxwell’s subsystem, care 
must be exercised to ensure that all the terms appearing in the fluxes 

®em>F em ,G em are expressed as functions of the primitive variables 

Qem = {E„ E y , E„ B x , B y , B z , q 0 }’ (59) 

For illustration, the flux vector E em can be extracted from equation (47) as 

0 ' 

H,/e 0 
-H y /e 0 

0 [ (60) 
-E z 

E y 

Jx ; 

For fluids with linear polarization and magnetization, H z and H y are the same as 
in equations (46), while J x is given in equation (55). The flux vector Jacobian 
matrix A_..is obtained as 

0 0 0 0 0 0 0 

a 21 a 22 ® a 24 a 25 a 26 0 

a 31 0 a 33 a 34 a 35 a 36 0 

0 0 0 0 0 0 0 

0 0 -1 0 0 0 0 

0 1 0 0 0 0 0 

a 71 a 72 a 73 a 74 a 75 a 76 V x 

where the coefficients are 
a 2 i=-X E v y a 3i=-X E v z (62) 

a 22=X Ey x a 33=X E V x 
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Matrices B £M and C £M may be obtained in the same fashion as equation (61). 

After tedious algebraic manipulations [26], the vector of eigenvalues of the 
flux vector Jacobian matrix A £M is found as 
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(72) 


This means that the eigenvalues X 1 = X 4 = 0, while X 7 = v x . The remaining 
four eigenvalues can be obtained from the fourth order algebraic equation 

X 4 +a EM X 3 + Vem X 2 + Yem ^ + Sem = 0 (73) 

where the coefficients in the fourth order characteristic polynomial are 
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For illustrative purposes, the following are the eigenvalues in the case of one- 
dimensional EMHD flow where v y = v z = 0 and a 22 = a 33 and a 25 = 0. Hence 
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Since //— equals the speed of light in vacuum, it seems that for most 

practical applications the incoming and the outgoing electromagnetic waves will 
not be influenced by the fluid except in the situations where the fluid is very 
highly ionized or when the fluid moves with a speed comparable to the speed of 
light. In the case of a pure electro-magnetics without any fluid motion, 
polarization, magnetization, or electric charges (v = P = M = q 0 =0), these 

eigenvalues reduce to the eigenvalues of Maxwell’s equations for electro- 
magnetic fields in vacuum [35] 


X= 0, 
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After introducing the similarity transformation matrix S £M of the flux vector 
Jacobian matrix A EM , the eigenmatrix X EM corresponding to A EM becomes 

=EM ^E» ^E» ^B» ^B> V x] (92) 

where X E , X E , Xg, Xg are given by equations (78-81). 

For locally one-dimensional problems, wave propagation direction is well 
defined. For multi-dimensional problems, there is no unique direction of 
propagation, because the flux vector Jacobian matrices A EM ,8^,0^, cannot 



be simultaneously diagonalized. Therefore, characteristic boundary condition 
analysis allows that only one of these matrices (relating to only one coordinate 
direction) can be diagonalized at a time. 

In the case that the x-coordinate is in the main flow direction, premultiplying 
the equation (58) with the inverse of the similarity matrix, , gives 
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Here, vector H EM is given as 
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For the hyperbolic system, time dependent boundary conditions could be 
derived based on the principle that outgoing waves are described by 
characteristic equations, while the incoming waves may often be specified by a 
non-reflecting boundary condition [31,32,36]. Following this approach, the 
characteristic and non-reflecting boundary conditions at the inlet boundary x = a 
and at the outlet boundary x = b can be given by the i-th equation of the system 
(93) 
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Here, the left eigenvector S. is the i-th row of S^ M and 
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7.2 Characteristic-based boundary conditions for Navier-Stokes subsystem 
Similar derivations can be used to determine analytical expressions for the 
eigenvalues and the non-reflecting boundary conditions of the Navier-Stokes 
subsystem of the unified EMHD as shown by Dulikravich and Jing [26]. 




— 2 

Terms related to d,d andVT will not be considered in the evaluation of 
coefficients of the flux vector Jacobian matrix A NS since they are associated 
with first derivatives of velocity, v, or temperature, T. The flux vector Jacobian 
matrix A NS = 3E NS /3Q NS then becomes 
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The coefficients in this matrix are given in detail by Dulikravich and Jing [26]. 
Eigenvalue vector of the flux vector Jacobian matrix A NS is 
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which can be written as a diagonal eigenvalue matrix 


X =diag[v x , X*, \+, XJ, X+] 


(103) 


The eigenvalues X*,Xy> X+ , X* are obtained analytically by solving a fourth 
order characteristic polynomial (similar to equation 73) where 

a NS = a 22 + a 33 + a 44 (104) 
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so that the four eigenvalues are 
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with the coefficients given by equations of the type similar to equations (82-88). 

Characteristic waves defined by the Navier-Stokes equations in the EMHD 
system have a great dependency on both fluid dynamics and electro-magneto- 
dynamics, in particular, the electro-magnetic properties of the media and 
electro-magnetic field quantities. When electric and magnetic fields are absent, 
these eigenvalues reduce to the well-known eigenvalues of a classical Navier- 
Stokes system for Newtonian, incompressible flows. These eigenvalues are 
{ v x, v x, v x, v x +c,v x -c). Here, the equivalent local speed of sound is defined 

as c = Vv^+(P/p). 

Following Thompson’s approach [30,31], non-reflecting boundary conditions 
for the Navier-Stokes subsystem are hence formulated as follows. The 
characteristic form of Navier-Stokes subsystem influenced by the electro- 
magnetic effects is possible to write as 
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where the i-th equation is 
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and the new source vector is 
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Here, the left eigenvector SJ 1 ^ is the i-th row of S 
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Practical implementation of Thompson-type [31-33,36,26] non-reflecting 
boundary conditions deserves further comments. The essence of his approach is 
that one-dimensional characteristic analysis can be performed by considering 
the transverse terms as a constant source term. In order to provide well-posed 
non-reflecting boundary conditions in multi-dimensional cases, substantial 
modifications may be required to take into account the transverse terms at the 
boundaries [37,38]. It should be emphasized that physically there are cases 
where flow information propagates back from the outside of the domain into the 
inside through the boundaries by the incoming waves [39]. This fact makes it 
possible that building a perfectly non-reflecting (absorbing) boundary condition 
[40] might lead to an ill-posed problem. Under these circumstances, corrections 
may be needed to make them partially non-reflecting. 

7.3 Numerical integration of EMHD system 

It is often highly desirable to have a time-accurate unsteady solution to the 
governing EMHD equations. One numerical integration algorithm that could be 
used is an advanced form of the dual time-stepping technique, also called an 
iterative-implicit technique, originally developed by Jameson [41]. 

To create an instantaneous picture of the solution of the entire EMHD 
system at a given physical time, equation (43) must be driven to zero in its 
entirety, not, as is commonly done in time-marching techniques by driving only 
the physical time-dependent term to zero. To this end, a pseudo-time derivative 
is added to the EMHD system (equation 43) which can be rewritten as 



or as 
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where R is a composite of the spatial and source terms and is called the 
residual. Thus, given a physical time step the governing equations are time 
marched in pseudo time, x . Upon convergence, the right-hand side of equation 
(118) becomes zero and the solution at the desired physical time level, t, is 

obtained. Note that the pseudo-time dependent variable vector, Q, does not 
have to be the same as the physical time dependent variable vector, Q . 

An additional concern of great importance is that the system of equations 
develops zero terms in the pseudo-time dependent variable vector, Q, for 
incompressible fluids, fluids without electric charges, or systems in which the 
electric and magnetic fields are non-interacting. This poses significant 

problems for time-marching numerical solutions. This problem may be 

alleviated, however, by proper selection of pseudo-time dependent variable 
vector, Q, and through the use of matrix preconditioning. 

By premultiplying Q with a properly selected matrix, it is possible to 
directly control the system eigenvalues. This prevents development of zeros in 
the pseudo-time dependent variable vector, Q, and vastly improves iterative 
convergence rates over a wide variety of flow regimes (low and high Mach and 
Reynolds number combinations). The preconditioning matrix, U(Q) , for the 
EMHD system could be based on one developed by Merkle and Choi [42] for 
the Navier-Stokes system. The preconditioned EMHD system may be written as 


r'3Q = £ 3Q 

™ 3x — 3t 


(119) 


Equation (119) can be transformed to a body-conforming non-orthogonal 
curvilinear time-dependent (£, T|, t) coordinate system. A high order of 
accuracy is desired to properly resolve unsteady motions. A finite difference 
scheme using fourth order accurate spatial differencing and second order 
accurate physical time differencing could be used while the solution is advanced 
in pseudo-time using a four-stage Runge-Kutta scheme which is second order 
accurate for non-linear problems. Fourth order accuracy should be selected for 



the spatial derivatives based on extensive research completed by Carpenter et al. 
[30] which found that a Runge-Kutta advanced fourth order accurate scheme 
provided the best convergence and stability of higher order schemes at 
reasonable computational cost. Second order accurate differencing in physical 
time could be selected based on stability and convergence studies performed by 
Melson et al. [43] who found that for a Runge-Kutta advanced dual time- 
stepping scheme second order backward differencing provided the most stable 
physical time discretization while providing excellent resolution. The new 
physical time step could be treated implicitly in pseudo-time, while all old 
physical time steps and spatial derivatives could be treated explicitly. This is 
unlike Jameson’s early method [41] that treats both the physical time and the 
spatial derivative explicitly and causes a restriction on the maximum physical 
time step allowed. The discretized preconditioned system may be written as 

Q°=Q n (120) 
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Q n+1 =Q 4 (122) 

where m=l,2,3,... represents the physical time step, n=l,2,3,... represents the 
pseudo-time step, and i=l,2,3,4 is the Runge-Kutta stage number. Also, 

£ = 3Q/3Q and a, are the Runge-Kutta coefficients. Note that the physical 

time-dependent term on the right hand side of equation (120) is held constant 
for all four Runge-Kutta stages. 


8. SUBMODELS OF EMHD 

Until now, the numerical solutions of the unsteady three-dimensional EMHD 
flows that have been reported in the open literature [34-36] did not account for 
polarization or magnetization effects and did not involve charge density 
transport equation. The reason is that the complete unified EMHD system is 
very large having extremely complicated source terms and two extremely 



different time scales for the electro-magnetic fields and the flow-field. 
Consequently, a number of simplified versions of the EMHD system have been 
traditionally used in practical applications. These submodels can be grouped in 
two general categories: EHD models and MHD models [1 1-13,44]. 

From the unified EMHD model, it can be seen that the electromagnetic field 
is not the only cause of electric current and that the temperature gradient is not 
the only source of heat conduction as is commonly assumed. The electric field, 
magnetic field, and heat conduction may couple to produce charge motion and 
heat transfer. These couplings are called phenomenological cross effects and 
may be placed in four general categories: 1) thermoelectric, 2) 

galvanomagnetic, 3) thermomagnetic, and 4) second order effects [9, p. 161- 
163]. These categories are based on the source of the effect and each will be 
described in turn, as will be a comparison between classical EHD and MHD 
models and the unified EMHD theory. The comparison concentrates on 
similarities and differences between electro-magnetic force and electric current 
and heat conduction terms in the EHD, MHD, and EMHD models. The 
inadequacies of simple superpositioning of classical simplified models to fully 
describe the unified EMHD flows are also noted. 

Couplings between the temperature gradient and the electric field cause 
thermoelectric effects so that a temperature gradient in the material produces an 
electric current (Thompson effect), while applied electric field produces heat 
conduction in the material (Peltier effect). These two effects together are 
known as the Seebeck effect and form the basis for thermocouples. Also note 
that the a, term in the electric conduction current (equation 22) and the k 4 term 
in the heat conduction (equation 23) are the ohmic charge conduction and 
Fourier heat transfer, respectively. 

When the electric and magnetic fields are simultaneously applied but are not 
parallel, electric current (Hall effect) and heat conduction (Ettingshausen effect) 
perpendicular to the plane containing the electric and the magnetic fields are 
induced in the media. These effects are termed galvanomagnetic [9, p.161-163]. 

When the temperature gradient and the magnetic field are simultaneously 
applied but are not parallel, electric current (Nemst effect) and heat conduction 
(Righi-LeDuc effect) perpendicular to the plane containing the temperature 
gradient and the magnetic field are induced in the material. These effects are 
termed thermomagnetic. 

It should be noticed (equation 22) that the interaction of the average rate of 
deformation tensor and the electric field can also create the electric current, 
while the interaction of the material deformation tensor and the electric field can 
create the temperature gradient (equation 23). These piezo-electric and piezo- 
magnetic effects can further be enhanced if the material is non-isotropic. 


8.1 Classical electro-hydrodynamics (EHD) 

As mentioned previously, EHD flows are those in which magnetic effects 
may be neglected and charged particles are present, while only a quasi-static 
electric field is applied so that the magnetic field, both applied and induced, 
may be neglected [11]. One of the implied assumptions is that the flows are at 
non-relativistic speeds, although in astrophysical flows this assumption cannot 
be made [1]. Atten and Moreau [44] present a detailed coverage of classical 
EHD modeling and discuss the relative importance of terms in the force and 
electric current through stability analysis. With these assumptions, the 
Maxwell’s system reduces to [1 1] 

V-D = V-(eE) = q 0 ( 123 ) 

^ + V-I=0 (124) 


With classical EHD assumptions, the electro-magnetic force in the unified 
EMHD theory reduces to: 

=q 0 E + (VE)P = q 0 E + (VE)-8 p E (125) 


This is not the form of the electro-magnetic force usually seen in classical EHD 
formulations [11]. Through the use of thermodynamics and the material 
constitutive equation of state, the electric force per unit volume in EHD is most 
often used in the following equivalent forms [10, p.505-507] [8, p.59-63] 
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(127) 


The three terms in the equation are the electrophoretic, dielectrophoretic and 
electrostrictive terms, respectively. 

The electrophoretic force or Coulomb force is caused by the electric field 
acting on free charges in the fluid. It is an irrotational force except when charge 
gradients are present [45]. 



The dielectrophoretic force is also a translational force, but is caused by 
polarization of the fluid and particles in the fluid. The dielectrophoretic force 
will occur where high gradients of electric permittivity are present. This 
condition will be true in high temperature gradient flows, multi-constituent 
flows, particulate flows [18] or any time the electric field must pass through two 
contacting media of different permittivities [46]. Grassi and DiMarco [47] treat 
the dielectrophoretic force as it applies to bubbly flows and heat transfer. 
Poulter and Allen [45] note that the dielectrophoretic force produces greatest 
circulation when the dielectric permittivity is inhomogeneous and non-parallel 
with the applied electric field. 

The last force, the electrostrictive force, is a distortive force (as opposed to 
the previous translational forces) associated with fluid compression and shear. 
The electrostrictive force is usually smaller than the -phoretic forces. It is 
present in high pressure gradient flows, compressible flows, and flows with a 
non-uniform applied electric field. Pohl [18] describes this phenomenon in 
greater detail. 

Classical EHD modeling derives directly from the unified EMHD theory. 
Thus, the electric current density, using EHD assumptions, reduces to 


I = q 0 v + a,E + a 4 VT (128) 

However, this is not the form seen in classical EHD models [11] which typically 
define the conduction electric current as only the first term of equation (21). 
However, more advanced classical EHD models define the current as [9, p.562] 

I = q 0 v + I c =q 0 v + q 0 bE-D 0 Vq 0 (129) 

The last two equations imply that the temperature gradient is directly related to 
the electric charge gradient. This may be shown to be true based on the 
Einstein-Fokker relationships, derived from studies of Brownian motion [25, 
p.264-273], which relate any concentration gradient to a charge mobility and a 
diffusion. Newman [48] also provides a detailed discussion of the concepts of 
diffusion and mobility. The electric charge diffusion term is often neglected 
where only limited amount of free charges are available [49]. 

By introducing classical EHD assumptions in the unified EMHD theory, the 
equation (23) for heat flux reduces to 


q = k,E + k 4 VT 


(130) 



The classical EHD models neglect the contribution to heat transfer from the 
electric field so that equation (130) reduces to Fourier’s law of heat conduction. 

q = -KVT (131) 

Although classical EHD modeling seems to neglects heat transfer induced by 
the electric field and electric current. Joule heating effect (- J c E term from 
EMHD equation 19) is usually included in the EHD computations [50,51]. 

8.2 Classical magneto-hydrodynamics (MHD) 

The classical modeling of MHD assumes non-relativistic and quasi- 
magnetostatic conditions. It implies that electric current comes primarily from 
conductive means and that there are no free electric charges in the fluid [11]. 
With these assumptions Maxwell’s system becomes 


VB = 0 

(132) 
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The modifications to the Navier-Stokes relations come from the electro- 
magnetic force on the fluid from which all induced electric field terms have 
been neglected. Using the MHD assumptions, the electro-magnetic force per 
unit volume in the unified EMHD theory becomes [11] 

£ em =IxB + (VB)M (136) 

The second term, source of dimagnetophoretic and magnetostrictive forces, is 
typically neglected in classical MHD [10, p.508]. Thus, the electro-magnetic 
force per unit volume in the classical MHD is modeled as [ 1 1 ] 

f EM = I x B 


(137) 



By making MHD assumptions, the conduction current in the EMHD can be 
expressed with equation (38). However, classical MHD theory usually defines 
the electric conduction current as [10, p.5 10] 

Ic =cy 1 E + a 4 VT = a 1 E + a,(vxB) + a 4 VT (138) 

Here, a 4 is the Seebeck coefficient [9, p. 174] which in some classical MHD 

formulations is not used [11]. Clearly, the classical MHD formulations neglect 
a significant number of physical effects [52,53]. 

Similarly, in classical MHD modeling, Joule heating is often included in the 
energy relation, but the heat transfer constitutive relation remains the same as in 
equation (11). In comparison, the unified EMHD model with classical MHD 
assumptions can be expressed with equation (39). 

It could be concluded that classical EHD models include many important 
effects and correspond to the unified EMHD theory well, while classical MHD 
formulations need improvements in the force, current and heat transfer terms. 

As in classical EHD modeling, it is important to be aware of the fact that 
many force, current and heat transfer terms can be written in several different 
forms, each of which is equivalent. It is, therefore, important to recognize the 
potential danger of simply adding terms from different MHD models. 


9. SOLIDIFICATION WITH ELECTRO-MAGNETIC FIELDS 

During solidification from a melt, if the control of melt motion is performed 
exclusively via an externally applied variable temperature field, it will take 
quite a long time for the thermal front to propagate throughout the melt thus 
eventually causing local melt density variations and altering the thermal 
buoyancy forces. It has been well known that an externally applied steady 
magnetic or electric field can, practically instantaneously, influence the flow- 
field vorticity and change the flow pattern in an electrically conducting fluid 
[51-59,33]. Similarly, it is well-known that applying an electric potential 
difference to a flow-field of a homogeneous mixture will cause fractionation or 
separation of the homogeneous mixture into regions having high concentration 
of the constituents. This phenomena, known as free-flow electrophoresis, has 
been extensively studied experimentally and, to a lesser extent, numerically [50] 
using classical EHD modeling. Nevertheless, there are no publications yet on 
actual algorithms for determining the proper variation of intensity and 
orientation of the externally applied magnetic and electric fields. This is not a 



trivial problem because we are dealing with a moving electrically conducting 
fluid within which an electric current is induced as the fluid cuts through the 
externally applied magnetic field lines [11]. This induced electric current 
generates heat (Joule effect) as it passes through the fluid that has a finite 
electrical resistivity. In the case of solidification, the amount of heat generated 
through the Joule effect due to the externally applied magnetic field is often 
neglected compared to the latent heat of solidification and the amount of heat 
transferred in the melt by thermal conduction. 

The latent heat released or absorbed per unit mass of mushy region (where 
Tiiquidus > T > T s0 [ idus ) is proportional to the local volumetric liquid/(liquid + 
solid) ratio often modeled [59] as 
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(139) 


Here, 0 is the non-dimensional temperature, the exponent n is typically 0.2 < n 
< 5, subscripts l and s designate liquid and solid phases, respectively, while f = 
1 for T>T liquidus and f = 0 for T<T solidus . Physical properties (density, 

viscosity, heat conductivity, heat capacity, etc.) are often significantly different 
in the melt as compared to the solid phase. We will assume linear variation of 
density as a function of the non-dimensional temperature, 0 , in the liquid 
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with a similar expression for the solid phase where the reference values are 
designated with the subscript "r". In this work, we assumed that electric 
conductivity and magnetic permeability do not vary with temperature. 

The EHD and the MHD systems of equations including solidification can be 
non-dimensionalized in a number of ways. The typical non-dimensional 
numbers are [33,60]: 
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where ji vr ,c r , A<j> r ,K r ,p r ,L r ,^ r are the reference values of viscosity, specific 
heat, electric potential difference, heat conductivity, magnetic permeability, 
latent heat of liquid-solid phase change, and length, respectively. Also, mixture 
density and modified heat capacity can be defined as 
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An enthalpy method [58,59] can be used to formulate the equivalent specific 
heat coefficient in the solid phase defined as 



J_3L 
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(147) 


so that latent heat is released in the mushy region according to equation (139). 


9.1 EHD and solidification 

EHD equations for phase-changing liquid-solid mixtures, where the solid 
phase is treated as the second liquid with extremely high viscosity, can be 
derived using Boussinesq approximation for thermal buoyancy [61]. We can 
also define mixture electric charge mobility 



( 148 ) 


bmix =f b # +(l-f)b s 
and combined hydrodynamic and hydrostatic pressures in liquid and solid 
p , 9 , p cp 

P'=— P s =- + 73- (149) 

P< Pr Ps Fr 

where cp is the non-dimensional gravity potential defined as g = Vcp . 

Assuming equal velocities for both phases, the mass conservation is 

V I = 0 (150) 

Linear momentum conservation for two-phase EHD flows with thermal 
buoyancy and Coulomb force 
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Energy conservation for incompressible two-phase EHD flows including Joule 
heating can be written as [60] 
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Electric charge conservation equation including migration and diffusion is 
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Since E = -V<(), the electric potential equation resulting from equation (13) 


V-[(fe, + (l-f> s )V<|>]=-N E q 0 


(154) 


must be solved simultaneously with the equations (151-154). 

9.2 MHD and solidification 

MHD two-phase solid-liquid flows can be modeled using a similar approach. 
The non-dimensional Navier-Stokes equations for phase-changing mixtures of 
two liquids (solid phase is treated as the second liquid with extremely high 
viscosity), can be formulated [33] so that the mixture mass conservation is 

Vv = 0 (155) 

Linear momentum conservation for two-phase MHD flows with thermal 
buoyancy and magnetic force 
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The non-dimensional hydrodynamic, hydrostatic, and magnetic pressures were 
combined to give 
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where cp is the non-dimensional gravity potential defined as g = Vcp. Then, the 

energy conservation for incompressible two-phase MHD flows including Joule 
heating can be written as [33] 
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The magnetic field transport equation for the two-phase MHD flow in its non- 
dimensional form becomes [1, p. 150] 
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If electric conductivity and magnetic permeability are assumed constant, then 
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needs to be solved either simultaneously or intermittently [33] with the 
equations (155-158). 
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